%% =======================================================================
% TABLE 4
% ========================================================================
clear all; close all; clc;

currentFile = mfilename( 'fullpath' );
[pathstr,~,~] = fileparts( currentFile );
addpath( fullfile( pathstr, '_aux' ) );
addpath( fullfile( pathstr, '_tables' ) );

rng(11);
NN    = 1e5;

mu        =  0;
tt_save   =  1/15;
tx_save   = [0.001:0.005:0.05];
tx1_save  = [0.01:0.05:1/2];
teps_save = [1.00:0.10:2];

data_moments = [-0.19,-0.20,0.41]; 

b_save    = zeros(length(tt_save),length(tx_save), length(tx1_save), length(teps_save));
beta_save = zeros(length(tt_save),length(tx_save), length(tx1_save), length(teps_save));
delta_save= zeros(length(tt_save),length(tx_save), length(tx1_save), length(teps_save));

for ii=1:length(tt_save)
    tt = tt_save(ii);
    for jj=1:length(tx_save)
        tx = tx_save(jj);

        for mm=1:length(tx1_save)
            tx1 = tx1_save(mm);
            
            if tx1<tx
                 beta_save(ii,jj,:,:)  = zeros(length(tx1_save),length(teps_save)); 
                 delta_save(ii,jj,:,:) = zeros(length(tx1_save),length(teps_save));
                 b_save(ii,jj,:,:)     = zeros(length(tx1_save),length(teps_save));
            else
            for nn=1:length(teps_save)
                teps = teps_save(nn);
                
        
        theta = mu*ones(NN,1)+sqrt(1/tt)*randn(NN,1);

        v     = tx1/(tt+tx1);
        vs    = tx/(tt+tx);
        ty    = vs^2;

        ei1   = sqrt(1/tx)*randn(NN,1);
        xi1   = theta+ei1;

        wx    = tx1/(tt+2*tx1+ty*teps);
        wy    = vs*teps/(tt+2*tx1+ty*teps);

        ei2   = sqrt(1/tx)*randn(NN,1);
        eps   = sqrt(1/teps)*randn(NN,1);
        xi2   = theta+ei2;
        y     = v.*theta+eps;

        fi1  = v.*xi1+(1-v).*mu;
        fi2  = wx.*xi1+wx.*xi2+wy.*y+(1-2*wx-wy).*mu;

        error = theta-fi2;
        rev   = fi2-fi1; 

        beta_save(ii,jj,mm,nn)  = rev\error; 
        delta_save(ii,jj,mm,nn) = y\error;   

        f1   = v.*theta+(1-v).*mu;
        f2   = wx.*theta+wx.*theta+wy.*y+(1-2*wx-wy).*mu;

        dev1  = fi2-f2;
        dev0  = fi1-f1;

        b_save(ii,jj,mm,nn)   = dev0\dev1;
        
            end
        
            end
        end
    end
end

%% Match
beta_save_net  = abs((beta_save - data_moments(1))./data_moments(1));
delta_save_net = abs((delta_save - data_moments(2))./data_moments(2));
b_save_net     = abs((b_save - data_moments(3))./data_moments(3));

WW = beta_save_net + delta_save_net + b_save_net;

[v,loc] = min(WW(:));
[II,JJ,MM,NN] = ind2sub(size(WW),loc);

disp('Regression Estimates'); 
disp([b_save(II,JJ,MM,NN) beta_save(II,JJ,MM,NN) delta_save(II,JJ,MM,NN)]);

NNN    = 1e6;

tt    = 1.00; 
tx    = tx_save(JJ)/tt_save(II);
tx1   = tx1_save(MM)/tt_save(II);
teps  = teps_save(NN)/tt_save(II);
mu    = 0;
theta = sqrt(1/tt)*randn(NNN,1);

v     = tx1/(tt+tx1);
vs    = tx/(tt+tx);
ty    = vs^2;

ei1   = sqrt(1/tx)*randn(NNN,1);
xi1   = theta+ei1;

wx    = tx1/(tt+2*tx1+ty*teps);
wy    = ty*teps/(tt+2*tx1+ty*teps);

ei2   = sqrt(1/tx)*randn(NNN,1);
eps   = sqrt(1/teps)*randn(NNN,1);
xi2   = theta+ei2;
y     = (v/vs).*theta+(1/vs).*eps;

fi1   = v.*xi1+(1-v).*mu;
f1bar = v.*theta + (1-v).*mu;
fi2   = wx.*xi1+wx.*xi2+wy.*y+(1-2*wx-wy).*mu;

err   = theta-fi2;
dev   = fi2-f1bar;
clem  = dev\err; 

%% Table
names        = {'Data', 'Model'};
cols1        = [-0.19; beta_save(II,JJ,MM,NN)];
cols2        = [-0.20; delta_save(II,JJ,MM,NN)];
cols3        = [0.41; b_save(II,JJ,MM,NN)];
cols4        = [-0.59; clem];
cols5        = [nan(1); sqrt(tx_save(JJ)/tt_save(II))];
cols6        = [nan(1); sqrt(teps_save(NN)/tt_save(II))];
cols7        = [nan(1); sqrt(tt_save(II)/tt_save(II))];
cols8        = [nan(1); sqrt(tx1_save(MM)/tt_save(II))];

table_new             = table(names',cols1, cols2, cols3, cols4, cols5, cols6, cols7, cols8); 
writetable(table_new,'_tables/table_e1.txt','Delimiter',' ')
table_new